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ABSTRACT 

We present the extension of previous two-dimensional simulations of the time- 
dependent evolution of non-relativistic outflows from the surface of Keplerian 
accretion disks, to three dimensions. As in the previous work, we investigate the 
outflow that arises from a magnetised accretion disk, that is initially in hydro- 
static balance with its surrounding cold corona. The accretion disk itself is taken 
to provide a set of fixed boundary conditions for the problem. 
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We find that the mechanism of jet acceleration is identical to what was estab- 
lished from the previous 2-D simulations. The 3-D results are consistent with the 
theory of steady, axisymmetric, centrifugally driven disk winds up to the Alfven 
surface of the outflow. Beyond the Alfven surface however, the jet in 3-D becomes 
unstable to non-axisymmetric, Kelvin-Helmholtz instabilities. The most impor- 
tant result of our work is that while the jet is unstable at super- Alfvenic speeds, 
it survives the onset of unstable modes that appear in this physical regime. 

We show that jets maintain their long-term stability through a self-limiting 
process wherein the average Alfvenic Mach number within the jet is maintained 
to order unity. This is accomplished in at least two ways. First, poloidal magnetic 
field is concentrated along the central axis of the jet forming a "backbone" in 
which the Alfven speed is sufficiently high to reduce the average jet Alfvenic Mach 
number to unity. Second, the onset of higher order Kelvin-Helmholtz "flute" 
modes (m > 2) reduce the efficiency with which the jet material is accelerated, 
and transfer kinetic energy of the outflow into the stretched, poloidal field lines 
of the distorted jet. This too has the effect of increasing the Alfven speed, and 
thus reducing the Alfvenic Mach number. The jet is able to survive the onset of 
the more destructive m = 1 mode in this way. 

Our simulations also show that jets can acquire corkscrew, or wobbling types 
of geometries in this relatively stable end-state, depending on the nature of the 
perturbations upon them. Finally, we suggest that jets go into alternating periods 
of low and high activity as the disappearance of unstable modes in the sub- 
Alfvenic regime enables another cycle of acceleration to super-Alfvenic speeds. 

Subject headings: Winds: accretion, accretion disks- proto-stars: jets-ISM: jets 
and outflows-MHD: magnetic fields 

1. INTRODUCTION 

Astrophysical jets are observed to be associated with young stellar objects (e.g. reviews 
by Konigl & Pudritz, 2000), as well as stellar-mass (Mirabel & Rodriquez 1998) and super- 
massive black holes (Begelman, Blandford & Rees 1984). Most, if not all, of the jets in this 
diverse collection are observed to have accretion disks associated with their central objects. 
Observations of jets from young stellar objects (YSOs) reveal that the accretion rate through 
the underlying disk is proportional to the mass loss rate carried in the jet, suggestive of a 
direct physical link between them. Moreover, both the radiation fields and rotation rates of 
low mass YSOs are observed to be far too feeble to launch either radiatively, or magneto- 
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hydrodynamically (MHD) driven winds, again suggesting that the outflow originates on the 
disk. The MHD nature of jets in YSOs was supported by the discovery of strong magnetic 
fields associated with the jet from T-Tauri (Ray et al. 1997), as well as by the observed 
narrowing of the opening angle of such jets with increasing distance from the central object 
in spite of dwindling external gas pressure (Burrows et al. 1996). 

Blandford & Payne (1982, henceforth BP) were the first to show that accretion disks, if 
threaded by large-scale, open magnetic field lines, will have their surface layers stripped away 
and flung out by the centrifugal stresses that act along such field lines. These outflows are 
subsequently collimated by the inescapable hoop-stress that arises from the toroidal magnetic 
field that develops within the jet itself. Jets in this picture are ultimately powered by the 
gravitational binding energy released during disk accretion. This simple, elegant picture of 
centrifugally accelerated disk winds potentially provides a universal model for jets in quasars 
(BP), microquasars (e.g., Mirabel & Rodriguez 1998) as well as YSOs (Pudritz & Norman 
1983, 1986). Furthermore, it is unlikely that jets are merely exotic oddities in astrophysics. 
MHD jets probably play a fundamental role in the physics of star formation, as well as black 
hole building because they can be even more efficient in stripping out the angular momentum 
of gas in the accretion disk than MHD turbulence (e.g., Pelletier & Pudritz 1992). 

Although the mechanism for the acceleration and collimation of MHD disk winds is 
conceptually clear, a detailed mathematical understanding of this phenomenon has proven 
to be elusive. Not much of a general quantitative nature is known beyond the predictions 
of conservation laws that pertain to time-dependent, axisymmetric, ideal MHD flows, or the 
predictions of models with additional simplifying assumptions (e.g., the imposition of self- 
similarity for time-independent, axisymmetric outflows, as in BP; Ostriker 1997; Ferreira 
1997). In spite of the concerted efforts of a large number of theorists, the challenge of 
finding general solutions to the highly non-linear, "Grad-Shafranov" equations for MHD jets 
(which have three types of MHD wave propagation) has proven to be insurmountable (see 
Heinemann & Olbert 1978, for a general discussion). A more general and natural approach 
to finding both stationary as well as time-dependent solutions to the jet problem therefore, 
is through the use of time-dependent MHD simulations. The advent of MHD codes over the 
last decade is rapidly transforming our knowledge about this rich and complex problem. 

Most of the simulations of MHD disk winds published to date are axisymmetric. The 
results of such simulations are in broad agreement with theoretical work on centrifugally 
driven outflows (e.g., Ouyed & Pudritz 1997a; 1997b; 1999, hereafter OPI, OPII, and OPIII; 
Romanovaet al. 1997; Kudoh, Matsumoto, & Shibata 1998; Meier et al. 1997; Krasnopolsky, 
Li, & Blandford 2002, hereafter KLB). The collimation of jets by the "hoop stress" engendered 
by their toroidal fields has been observed in simulations that posit modest gradients of the 
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magnetic field across the underlying accretion disk (e.g., OPI; KLB). One caveat is that 
outflows produced by disks threaded by a steeply declining poloidal magnetic field, such as 
the extreme case of the split-monopole distribution of Romanova et al. (1997), do not appear 
to collimate well. Spruit, Foglizio & Stehle (1997) argue that jets from accretion disks should 
be strongly unstable to helical modes, and that jets, therefore, might be collimated by larger- 
scale poloidal magnetic fields rather than by magnetic hoop stress. 

The purpose of the present paper is to investigate one of the principle remaining chal- 
lenges in the theory of jets, namely, why are real 3-D jets so stable over great distances 
in spite of the fact that they are probably threaded by strong toroidal fields? It is well 
known that the purely toroidal field configurations that are used to help confine static, 3-D, 
Tokomak plasmas are unstable (e.g., Roberts 1967; Bateman 1980). The resulting kink, or 
helical (m — 1) mode instability derived from a 3-D linear stability analysis is powered by 
the free energy in the toroidal field, namely B^/8n (Eichler 1993). Shouldn't this instability 
affect the longevity of astrophysical jets too? 

There are several major differences between the physics of astrophysical jets, and the 
simpler Tokomak configurations. First, jet plasmas are not static but consist of gas that 
typically has been accelerated to velocities greater than the fast magnetosonic (FM) wave 
(e.g., with an FM Mach number, Mfm, of several). A jet moving with sufficiently large FM 
Mach number, and which carries a current, can also suppress the growth of several types 
of instabilities. Second, astrophysical jets are threaded by a "backbone" of purely poloidal 
magnetic flux which may act to stiffen the jet against short-wavelength kink instabilities. 
One anticipates therefore that jets with poloidal fields at their core might still be unstable 
to wavelengths much longer than their width, but may very well survive such transverse 
motions. Finally, in previous simulations of 2-D jets (OPI), it was found that in the final, 
stationary jet, twice as much energy is carried by the bulk poloidal outflow (kinetic energy) 
than by the dominant toroidal magnetic field. This may also be a factor that favours the 
stability of jets, if it can be demonstrated to occur in 3-D outflows as well (see also KLB). 

Analytic studies of the stability of 3-D, non-relativistic jets have focussed mainly on 
linear stability analysis of rather idealised jet configurations. These assume that jets have a 
radial, top-hat velocity profile, such as that of the pressure driven flow arising from an over- 
pressured orifice. The stability of 3-D jets of this type, which also contain purely poloidal 
or toroidal magnetic fields, is discussed by Ray (1981); Ferrari, Trussoni & Zaninetti (1981); 
Fiedler & Jones (1984), and many others. As an example, Ray (1981) showed that the 
growth rate of the Kelvin-Helmholtz (K-H) helical (kink) modes for all wavelengths longer 
than the jet radius R (i.e., kR < 1), vanish if flows are sub-Alfvenic. 

The growth rates of helical modes in super-Alfvenic jets is lower than that of purely 
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hydrodynamic jets. More recent work by Appl & Camezind (1992) examined the stability 
of more general, current-carrying jets that contain force-free, helical magnetic fields. These 
latter systems are more stable than their purely hydrodynamic, or even purely poloidally 
magnetised, counterparts. These authors also found that jets are increasingly stabilised at 
progressively larger scales as the jet Mfm is increased. 

Several extensive numerical simulations and studies of the stability of simple, 3-D, uni- 
form jets may be found in the literature. For example, Hardee & Rosen (1999) performed 
3-D simulations of "equilibrium" jets, and found that these uniform, magnetised jet models 
remain Kelvin-Helmoltz stable to low-order, surface helical and elliptical modes (m = 1,2), 
provided that jets are on average sub-Alfvenic. This is in accord with the prediction of 
linear stability analysis. However, most configurations for jet simulations use rather ad hoc 
prescriptions for the initial toroidal field configuration (e.g., Todo et al. 1993; Lucek & Bell 
1996; Hardee & Rosen 1999, to cite only a few) so that it is difficult to assess how pertinent 
the results are to the case of a jet that establishes its own toroidal field as the jet is accel- 
erated from the accretion disk. In general, the available analytic and numerical results for 
the stability of simple jets show that the fastest growing modes are of K-H type. These K-H 
instabilities are increasingly stabilised for super-Alfvenic jets, as Mfm is increased much 
beyond unity. It is also generally known that sub-Alvenic jets are stable. Taken together, 
these results suggest that 3-D jets are the most prone to K-H instabilities a bit beyond 
their Alfven surface 3 , a region wherein their destabilising super-Alfvenic character cannot 
yet be offset by the stabilising effects engendered at large super FM numbers (e.g., Hardee 
& Rosen 1999). If this is so, then simulations should include the region not too far from the 
acceleration zone, beyond the putative Alfven surface. 

This paper presents a numerical investigation of the stability of 3-D MHD jets launched 
by accretion disks. The strategy is to extend to three dimensions previous work on the simu- 
lation of axisymmetric, time- dependent jets from Keplerian accretion disks (OPI, OPII, and 
OPIII). These 2-D simulations made use of the ZEUS-2D code of Stone & Norman (1992a, 
1992b) and demonstrated the existence of stationary jets with properties that well match 
the predictions of the theoretical literature. They also revealed a class of new, intrinsically 
time-dependent episodic jets. It is natural therefore to extend our basic model approach, 
which rests on a secure numerical and physical foundation, to 3-D. Our basic finding is that 
3-D jets, driven by an underlying disk for a particular magnetic configuration, are ultimately 
stable. The simulations show that there are mechanisms that help preserve the integrity of 
real astrophysical jets despite the growth of unstable modes within them. 



3 The Alfven surface is denned by the set of Alfven points where Ma — 1, {ta, za), of the flow along each 
jet field line. 
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We outline the physical model underlying all of our simulations in §2 and detail our 
numerical methods in §3. We show representative simulations from two classes of simulations 
in which the Kepler rotation of the disk (a boundary condition of the simulations) is either 
sharply or gradually truncated at the inner disk edge (sections §4 and §5 respectively). We 
then analyse our data and show that stabilisation of jets by non-linear saturation of K-H 
instabilities appears to be the basic mechanism. We conclude in §6. 



2. BASIC PHYSICAL MODEL 



In an earlier series of papers (OPI, OPII, and OPIII), the behaviour of time-dependent, 
non-relativistic disk winds in 2-D was investigated with the intent of studying their accelera- 
tion and collimation. To this end, a particularly simple, and analytically tractable model for 
an accretion disk and its surrounding corona was chosen with a common threading magnetic 
field. The key simplification in this approach is to examine the physics of the outflow for 
fixed physical conditions in the accretion disc (see also Romanova et al. 1997; Ustyugova et 
al. 1998). In part, this simplification may be justified by the fact that typically, accretion 
discs will evolve on longer time scales than their associated jets. We retain this approach in 
the present 3-D work. 

Our physical model consists of an accretion disc whose surface pressure matches the 
pressure of an overlying, non-rotating corona in hydrostatic balance within an unsoftened 
gravitational potential well from a central object. As discussed in OPI, the resulting analyt- 
ical model for the corona in scaled units is then given by: 



where p is the matter density, $ is the gravitational potential, r is the radial distance from 
the central object, and 7 = 5/3 is the adiabatic index of the gas. 

In OPI and OPII, the pressure was given by assuming a strict polytropic gas, with 
polytropic index equal to 7 (thus, p ~ p 7 ), and we follow this strategy here. 4 Therefore, no 



4 This has the advantage of simplicity, particularly when trying to establish a numerically stable atmo- 
sphere using the gravitational potential. However, the disadvantages of using a strict polytrope include 
not being able to track contact discontinuities (e.g., between the disc and the corona) and not getting the 
Rankinc-Hugoniot jump- conditions right (entropy is strictly conserved everywhere, including across shocks). 
Independent simulations performed by D.A.C. using the full energy equation show that these concerns have 
minimal impact on the simulations presented herein because much of the flow is subsonic. Qualitative differ- 
ences between the polytropic and adiabatic equations of state appear only for supersonic flows where shocks 
and contact discontinuities play dynamically significant roles. 
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separate internal energy variable needs to be tracked. 

A corona in hydrostatic balance is implemented numerically so that left unperturbed, 
it remains in perfect numerical balance indefinitely to within machine round-off errors. We 
note in passing that without due care, the atmosphere can be extremely unstable numeri- 
cally, particularly when the internal energy equation is used. In particular, if numerical errors 
generated near the origin where the gradients are extremely steep are not specifically ad- 
dressed, interpolation errors will cause pressure gradients near the origin to be consistently 
overestimated, thereby launching a thin, fast jet along the symmetry axis. These jets are 
completely numerical in origin, destroy the atmosphere through which they propagate before 
the physical outflow can be established, and at the very least, corrupt the results from the 
physical jets launched later by fully physical means. Additional details will be published in 
a future paper. In part, our use of the polytropic equation of state was to help squelch this 
very effect. 

As discussed in OPI, the physical conditions on the rotation axis r = of this model 
require special attention. The accretion disk will have an inner edge, either because it abuts 
the outer surface of a YSO, or because it is terminated by the magnetosphere of the central 
object. Any gas or objects interior to this inner disk edge are taken to be non- rotating. 
In this paper, we accomplish this in two different ways. First, as in OPI and OPII, we 
simply truncate the Keplerian velocity profile abruptly at the inner radius (simulations A- 
D). Second, in simulation E, we adopt a velocity structure on the inner disk edge that falls 
to zero in a less dramatic way, as one might expect in the presence of a boundary layer (see 
Appendix A, §A.l). 

Some authors have posited the existence of an inner jet on the r = axis (e.g., KLB), 
rather than set up an initial hydrostatic state. This is done to mimic the possible existence 
of a jet from a rotating central object in its own right. As already discussed, we have seen 
that numerically driven jet-like outflow can arise if care is not exercised in setting up an 
initial hydrostatic state. Given the sensitivity of all calculations to boundary conditions on 
the axis, we feel it is imperative to establish the existence of an inner jet by a calculation 
that specifically includes all the complexities of a central, magnetised, rotating body, rather 
than just assuming it, or possibly creating it with numerical truncation errors. 

The external MHD torque on the disk, ultimately responsible for the disk wind, requires 
that a plausible threading magnetic field configuration be introduced. There are many 
possible choices, and we select among those that can be easily initialised within the ZEUS 
framework. In the previous 2-D work, current-free (and therefore force-free) configurations 
(J = 0), were used, so as not to disturb the equilibrium of the hydrostatically stable corona. 
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The most obvious choice is a constant vertical magnetic field (B z = constant), as used 
in OPII. OPI described the results of launching a jet into a non-uniform (but still force-free) 
magnetic field, consistent with what one would expect from a conducting plate in a vacuum 
in axisymmetry (Cao & Spruit 1994). Because the magnetic field lines in the configuration 
do not follow the grid axes, it is necessary first to evaluate the 0-component of the vector 
potential (A^), and then from this, determine the r- and ^-components of the magnetic field 
(technical details for a similar problem are given in Appendix A, §A.3). Only in this way 
can the initial magnetic field be established with a zero divergence everywhere on the grid 
to within machine round-off errors. 

Given that the 3-D simulations in the present paper are performed on a Cartesian, 
rather than a cylindrical co-ordinate system (see §3.2), we chose to simulate the behaviour 
of jets launched in an initially uniform magnetic field, parallel to the vertical z axis. This 
is far simpler than implementing the axisymmetric, potential configuration of Cao & Spruit 
(1994) in 3-D because of the difficulties in having to deal with three components of the vector 
potential 5 . Further, as seen in OPI and OPII, the vertical magnetic field model is, in many 
ways, more interesting. 

2.1. Parameters of the Model 

Our physical model treats the accretion disk as a set of fixed boundary conditions for 
the atmosphere. Thus, the values of five flow variables must be prescribed at every point of 
the accretion disk surface at all times (see also KLB). These five variables include the mass 
density, p(r); two components of the magnetic field B, namely B z (r) and B^r) [B r (r) is 
fixed by the solenoidal condition]; and finally, two components of the velocity field, namely 
v z (r) and v^r) [v r (r) in the disk is taken to be zero]. Note that the product pv z prescribes 
the mass loss rate per unit area from the disk, or, equivalently, the mass loading of field 
lines, and this quantity remains fixed throughout our simulations. We chose the value of the 
vertical speed, v z , to be far less than the sound speed of the disk and the expected, slow 
magnetosonic velocity. Thus the resulting disk wind, if it achieves steady state, is accelerated 
through all three critical points (see OPIII for detailed discussion). We retain exactly this 
approach in establishing the disk conditions in our full 3-D problem. 

Therefore, the length scale, rj, density of the corona, pi, and Keplerian velocity of the 
disc, vx,i, are all taken to be unity at the inner edge of the disc. Thus, the time scale is in 



5 More elaborate magnetic field configurations can be implemented/initialised using the JETSET tool 
developed in J0rgensen et al. (2001). 
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units of U = Ti/vK,i- The rest of the variables are set by using the functional forms given by 
equations (1), and by setting five basic parameters introduced in OPI. These are discussed 
below. 

First, while we assume that the disc and corona are in pressure balance everywhere 
along their common boundary there is a density jump (and thus a contact discontinuity) 
and this is given by the first parameter rji, the ratio of the disc density to corona density at 
the inner radius (r = 1). 

Second, while not critical to the ultimate launching of the jet, we do introduce a force- 
free toroidal magnetic field in the disc only. In part, this is done to reduce the gradient in 
Bfy across the disc-corona boundary, once the coronal toroidal field is established. To remain 
force-free, we choose the form: 

a, = (2) 

where is the second parameter, and is equal to the desired ratio B^j B z at r — 1. 

Third, the vertical velocity v z which provides the mass loading of the coronal magnetic 
field, is given by: 

V z = VinjVfr (3) 

where is the rotation velocity profile (mostly Keplerian) of the disc (Appendix A, §A.l) 
and Vi n j is the third parameter. 

Fourth, the ratio of Keplerian to thermal energy densities is given by: 

^ = —7- = -, 4 

Pi I pi Pi 

since vk,% = 1 and pi — 1. It is well known that for thermal atmospheres in hydrostatic 
balance, this ratio is given by 5; = 7/(7 — l) = 5/2fora7 = 5/3 gas. However, OPI argued 
from observational data that only a small fraction of the total pressure can be thermal 
in origin, and the rest must come from some other isotropic mechanism, such as Alfvenic 
turbulence. OPI considered the total pressure to originate from two terms, one thermal, 
the other turbulent, but since they modelled the Alfvenic turbulent pressure with the same 
7 = 5/3 polytrope used for the thermal component, the two components became physically 
indistinguishable. Thus, whether one thinks of the total pressure as being thermal plus 
Alfvenic, or all thermal, the numerical solutions are identical. Therefore, in the present 
work, we retain 5, for consistency with OPI and OPII, but note that it simply refers to the 
inverse of the portion of the pressure designated as thermal. In any case, the total pressure 
at Ti is equal to (7 — l)/7- 
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Finally, the fifth parameter is the plasma beta, given as usual by: 

Pi = |H> (5) 

where p t ,i is the initial thermal pressure at rj, and where the factor of Hq (or 4-7T, depending 
on one's units of choice) has been absorbed into the scaling for the magnetic field. Since 
p t ,i = we have: 

B z , = y/2/SiPi. (6) 

In the five 3-D simulations presented herein (see §3.3), the values for these five parameters 
are the same as in OPII, namely: 

(r^/i,,^,^, A) = (100.0,1.0,0.001,100.0,1.0). (7) 



3. 3-D SIMULATIONS 

3.1. Computational Details 

The 3-D simulations were computed with ZEUS-3D, a multi-dimensional, finite-differ- 
ence, Eulerian MHD code developed by D.A.C. While sharing a common lineage with the 
NCSA's public-domain code of the same name, our code uses different algorithms for solving 
the induction and momentum equations, namely, the Consistent Method of Characteristics 
(CMoC, see Clarke 1996 for details). 

Like all codes in the ZEUS family, our version of ZEUS-3D uses a staggered mesh to 
reduce the number of interpolations required at the zone faces. Interpolations are performed 
with the second-order accurate scheme first proposed by van Leer (1974). It uses an operator- 
split, time-centred algorithm for transporting the variables and applying the source terms. 

The CMoC algorithm sets our code apart from its predecessors. It uses a planar-split 
strategy, rather than the traditional directional splitting employed by earlier versions of 
the code (e.g., Stone & Norman 1992a; Hawley & Stone 1995). In addition, the magnetic 
induction, momentum transport, and transverse Lorentz acceleration are all tightly coupled 

— * 

by using the same interpolated values of B and v. It is a robust algorithm, and possesses no 
known numerical instabilities or problems in 3-D. 

Three relatively minor modifications were made to the code in order to perform the 
simulations discussed in the next sections. First, in order to "turn off" the internal energy 
equation and implement the polytropic equation of state, we simply replaced the pressure 
gradient terms in the momentum equation with gradients of the function p 7_1 + $, where p 
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is the updated density distribution, and $ is given by equation (1) at all time steps. Since 
the internal energy equation contributes to the dynamics only via the pressure gradient in 
the momentum equation, this is all that has to be done to affect this change. For example, 
other than for computational efficiency, there is no need to prevent the code from updating 
the internal energy (e); it is simply never used. 

Second, since the internal energy variable is ignored, it is necessary to modify the 
dependence of the CFL time step on the sound speed. Thus, instead of computing c s from e 
and p [i.e., c 2 s = 7(7 — l)e/p], we use c 2 s = 7P 7 " 1 , thereby introducing the polytrope explicitly. 

Third, a subroutine is required to initialise the corona (all flow variables, including the 
gravitational potential) and the disc boundary according to the discussion in the previous 
section. 

Other changes of a more technical nature were required to address boundary condition 
problems, special graphics, and other such things. As mentioned in §2.1, the hydrostatic 
atmosphere is very prone to numerical instabilities, particularly at the boundary, and we 
found these problems to be even more significant in 3-D. These details are relegated to 
Appendix A. 

3.2. Initialising the Simulations 

Contrary to what may be intuitive, it is inadvisable to perform a 3-D simulation such 
as this in cylindrical coordinates. For one thing, special treatment must be introduced for 
the "wedge zones" that abut the z-axis (no longer a symmetry axis in 3-D), and velocities 
that pass through the z-axis pose a very difficult numerical problem. Second, even with such 
technical details solved, plane waves are badly disrupted upon passing through the z-axis 
(John Hawley, private communication), and this provides an undesirable bias to what should 
be an unbiased three dimensional calculation 6 . 

Thus, we use Cartesian coordinates, (x,y,z), for these simulations. The disc is taken 
to lie along the x-y plane, and the disc axis corresponds to the z-axis (see Fig. 1). While 
Cartesian coordinates are the natural system to use to avoid any directional biases, it does 
introduce some of its own problems not encountered in the 2-D cylindrically symmetric 
simulations, and we discuss these further in Appendix A. 

Five separate 3-D simulations were performed for this work, and their details are sum- 



6 By unbiased, we mean that no axis should be preferred numerically in any way over another. 
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marised in Table 1. Simulations A, B, C, and E were all performed at the same resolution and 
same spatial extent, while simulation D was performed with a slightly larger spatial extent 
and with twice the radial resolution as the other four runs. In this paper, we concentrate 
primarily on simulations D and E. 

In units of the inside radius of the disc, the 'primary computational domain of runs 
A, B, C, and E have dimensions (—15 : 15, —15 : 15, : 60), and is divided into (60, 60, 120) 
uniform cubical zones. Thus, there are only two zones between the disc axis and the inner 
radius of the disc. This should be compared with the 2-D runs of OPI and OPII which, using 
cylindrical coordinates (z,r), was computed over a domain (0 : 80, : 20) and resolved with 
(500, 200) zones (and thus ten zones per r,). By comparison, therefore, these 3-D simulations 
have one fifth the resolution of the 2-D runs. 

The primary computational domain of run D has the same physical extend as the 2-D 
runs in OPI, namely (-20 : 20, -20 : 20, : 80), and is divided into (100, 100, 160) zones. 
Along the z-axis, the zoning is uniform, and thus has the same axial resolution (two zones 
per rj) as the lower resolution runs. However, along the x- and y-directions, 40 uniform zones 
resolve the range ± 5rj (giving a radial resolution in the vicinity of the disc axis of 4 zones 
per r^, while 60 "ratioed" zones are used to resolve the remainder of the range (—20 : —5 
and 5 : 20). At r = 15, the radial extent of the zones has grown to about 0.6rj (comparable 
to the lower resolution runs), and at r = 20, to about 0.9rj. Thus, resolution far away from 
the disc axis has been sacrificed to some extent in favour of higher resolution in the more 
important regions nearer the disc axis. 

Figure 1 shows a schematic of the setup of our grid for runs A, B, C, and E. The 
primary computational domains of all simulations are embedded in a greater computational 
domain which, for runs A, B, C, and E, is (—30 : 30, —30 : 30, : 120), while for run D, 
is (—40 : 40, —40 : 40, : 160). The portions of the greater domains that lie outside the 
primary domains are resolved with 10 to 20 severely ratioed zones, and are never used for 
analysis. They are merely there to act as a "buffer" between the primary computational 
domain and the imperfect outflow boundary conditions (see Appendix A, §A.4 for further 
discussion). 

We use inflow boundary conditions at the z = boundary (disc), even inside the inner 
disc edge. However, because the azimuthal velocity profile, [equation (A. 2)], goes to zero 
inside r = and because the inflow velocity, v z , is set to a fraction of v^, actual inflow is 
restricted to the portion of the boundary where ^ 0, namely the putative disc. In addition, 
because the Cartesian grid is rectangular, numerical problems arise near the corners of the 
grid if rotation of the fluid is permitted to persist there (Appendix A, §A.l). Thus, the 
Keplerian profile of the disc is reduced smoothly to zero between an "outer radius" (r G ), 
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which roughly corresponds to the radius of the smallest cylinder that can fully contain the 
primary computational domain, and r max , chosen to be greater than r Q but still less than 
the maximum extent of either the x- or y-axes (see equation A. 2 in Appendix A). Note that 
because r c lies outside the primary computational domain, the truncation of the Keplerian 
profile at r D has no visible effect within the region of analysis (namely the primary domain). 

As described, the corona and disc possess quadrantal symmetry, and runs A, B, and 
C were designed to determine the best way of breaking this. Run A was set up with the 
identical parameters and boundary conditions as the 2-D run in OPII, without any deliberate 
attempt to break the quadrantal symmetry. In principle, 2-D slices through this run should 
be very similar to lower resolution 2-D runs from OPII, and until t = 150, this was indeed 
the case. However, quadrantal symmetry is not the same as cylindrical symmetry, and their 
differences show up in the latter half of the simulation. Still, the jet remained centred about 
the disc axis and propagated at roughly the same velocity as observed in 2-D. 

In run B, the quadrantal symmetry was broken by offsetting the centre about which the 
velocity profile of the disc is computed relative to the centre of the gravitational potential 
(located at the grid origin). Initially, the offset (one to tens of percent of r i: it does not seem 
to matter) is along the x-axis and at t=10, is "jerked" to the same position along the y-axis. 
Meanwhile, in run C, the disc is effectively wobbled, rather than being jerked. Details of 
how the jerk and wobble are implemented are found in Appendix A (§A.2). 

Qualitatively, runs B and C are identical, and the quantitative differences are slight. 
However, both simulations are completely different from run A, in which nothing deliberate 
was done to break the symmetry. Thus, we conclude that it does not really matter how the 
quadrantal symmetry is broken, so long as it is broken. Run D is the same as run C but at a 
higher resolution and is discussed further in §4. No further discussion will be given for runs 
A, B, and C. 

Run D uses the same radial profiles for the velocity, magnetic field, and density in the 
disc as OPII. In particular, at the inner radius of the disc, the Keplerian velocity profile 
is suddenly truncated, leaving a cusp in at r = (Fig. 2, see also Appendix A, §A.l). 
Run E, therefore, was designed to test the dependence of the numerical results on how the 
Keplerian velocity profile is truncated at r — r^. In particular, in run E, we use the region 
between r = and r = 2r, to round off the Keplerian velocity profile smoothly to zero 
(Fig. 2). Thus, in run E, there is no cusp in the profile for v$ and the amount of angular 
momentum transferred from the disc to the corona is significantly less than that of run D. 
We find that this gives qualitative differences in the relative importance of the various modes 
of the K-H instabilities that are excited (see §5). 
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The greater computational domains of runs A, B, C, and E include about 900,000 zones 
and required about 9,000 MHD cycles and 50 hours of CPU time on an IBM Power 3+ 
processor to run to t — 300 (in units of Ti/vK,i)- For contrast, the greater computational 
domain of run D included nearly 2.9 million zones and required about 22,000 MHD cycles and 
18 days on the same machine to run to t — 400. In all cases, the simulation was terminated 
once the working surface of the jet left the greater computational domain, and/or when the 
effects of the imperfect outflow boundary conditions made themselves felt on the primary 
computational domain. 

Finally, the differences between runs C and D were surprisingly slight, affecting primarily 
the details of the profile of the jet. Thus, we justify using the lower resolution setup for run 
E, and performing a simulation that took two days rather than eighteen. 



Simulation 


primary 


greater 


ro/ri 


^max/^j 


symmetry 


inner v$ 




domain (r-j) 


domain (r-j) 






breaking 


profile 


A 


30 x 30 x 60 


60 x 60 x 120 


22 


28 


none 


truncated 


B 


30 x 30 x 60 


60 x 60 x 120 


22 


28 


jerked 


truncated 


C 


30 x 30 x 60 


60 x 60 x 120 


22 


28 


wobbled 


truncated 


D 


40 x 40 x 80 


80 x 80 x 160 


30 


38 


wobbled 


truncated 


E 


30 x 30 x 60 


60 x 60 x 120 


22 


28 


jerked 


smooth 



Table 1. Specifics of the five simulations performed in this work. The numerical 
resolution of the primary domains for simulations A, B, C, and E was 60 x 60 x 120 
zones, while for simulation D, 100 x 100 x 160 zones. The "truncated" and 
"smooth" profiles for v<f, are shown in Fig. 2. 



4. A "CORKSCREW" JET 

We first focus on simulation D, which other than dimensionality and resolution, was 
initialised in the same way as the simulation discussed in OPII. 

The quantities illustrated in false colour in Figs. 3a and 3b are respectively: 

J p(l) dl, J V • v(l) dl, 

where p{l) and V • v(l) are the density and velocity divergence at coordinate / integrated 
along the line of sight. Thus, Fig. 3a is a 2-D map of the column density in the primary 
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computational domain, while Fig. 3b indicates regions of compression (blue) and expansion 
(red) along the line of sight. In both cases, the z-axis of the data cube has been rotated 
by 20° out of the visual plane. The image is taken at t — 320, and is representative of the 
appearance of the jet from t = 210 through t = 400. The disc (not visible in this rendering 
because boundary values are not included in the line-of-sight integrations) is at the left hand 
side of the image, and outflow is from left to right. 

Figure 3a shows that the jet has settled into a quasi-steady state structure in the shape 
of a single helix, or a "corkscrew", and the bulk of this section is devoted to explaining 
how the jet reaches this configuration. On the primary computational domain, there are 
roughly 1.2 wavelengths of the helix, the wavelength itself depending on the complexities of 
the K-H instabilities ultimately responsible for the structure. The line-of-sight integration of 
the velocity divergence (Fig. 3b) is included because it nicely illustrates the most dynamic 
portions of the helical jet (two narrow blue ribbons on either side of the material jet visible 
in Fig. 3a indicating regions of greatest compression) twisting around a relatively stable, 
slightly expanding core (red). One must be careful in interpreting the velocity divergence 
as shocked regions, since the polytropic equation of state precludes the Rankin-Hugoniot 
jump-conditions. Thus, Fig. 3b is included for illustrative purposes only. 



4.1. The Nature of the Outflow 

Figure 4 shows a time-series of eight contour plots of density taken along the y-z plane, 
where the z-axis (horizontal) corresponds to the axis of the disc. In these slices, the +x-axis, 
located at y = z = 0, points into the page. The eight epochs chosen are t = 20, 60, 100, 160, 
200, 240, 320, and 400. In Fig. 4a, the vertical lines mark the location of the cross- sectional 
cuts shown in Figs. 7-10. Figures 5 and 6 are similar montages for the normal magnetic field 
(B^) and poloidal velocity vectors respectively. 

As in 2-D, a global Alfven wave is launched from the disc as the rotation twists the ini- 
tially uniform and vertical magnetic field in the corona. By t — 20, this wave has propagated 
a third of the way across the grid (Fig. 5a) and the collimated outflow has reached about 
z — 10 (Fig. 4a). Until t = 100, the jet behaves much like the 2-D jet reported in OPII as 
the non-axisymmetric modes have not grown enough to break the initial cylindrical symme- 
try. Thus, knots are generated in the same way they were in 2-D in response to the early 
appearance of the m = pinching mode of the K-H instability 7 (e.g., the four symmetrically 



7 We remind the reader that m = K-H mode pinches jets ("sausage" instability), but does not disrupt 
them. This is the only mode that can be excited in 2-D axisymmetry. In 3-D however, helical modes (m = 1) 



-16- 



positioned knots near the head of the advancing jet, each labelled with an 'H' in Fig. 4c). 

By t — 120 (not shown), the weak m = knots forming along the jet axis have all 
but merged back together, as the differences between the cylindrical and near-quadrantal 
symmetries start to appear. In particular, while only the m = mode can appear in a cylin- 
drically symmetry system, all modes that are multiples of four will appear in a quadrantally 
symmetric system, and the m = 4 mode begins to dominate the m = mode after t = 100. 
Figures 7 through 10 show x-y profiles at z = 30 (indicated by the vertical line in Fig. 4a) for 
the density, Alfvenic Mach number, velocity, and magnetic field (in the latter two, contours 
indicate the normal component, vectors indicate the poloidal component). The quadrantal 
distortion (m = 4 mode) is clearly evident in panels d, e, and f of Figs. 8 and 9. 

By t = 160, the effects of breaking the quadrantal symmetry are evident, as a clear 
sinusoidal distortion finally appears along the jet axis (Figs. 4d, 5d, and 6d). This is the 
beginning of the m — 1 mode of the K-H instability, and represents the greatest threat to 
the ultimate stability of the outflow. 

A brief recap of the various modes of the Kelvin-Helmholtz instability encountered so 
far is in order. First, Hardee, Clarke & Rosen (1997) show that while the m = mode is 
the fastest growing mode, its amplitude is always less than that of the m = 1 mode, whose 
amplitude is greater and whose growth rate is faster than all the other modes (m > 2). 
Thus, it is not surprising that we should first see hints of the m = mode discussed in OPII 
before the onset of the m = 1 mode. What may be surprising, however, is the appearance 
of the m = 4 mode before the m = 1 mode. 

In fact, the development of the two modes arises from very different processes. The 
m — 1 mode arises entirely from the growth of signals propagating from the left-hand 
boundary (e.g., the disc wobble). On the other hand, the m = 4 mode arises from in 
situ perturbations associated with the fact we are resolving a rotating, initially cylindrically 
symmetric object on a Cartesian grid. Thus, the m = 4 mode appears well before the m — 1 
mode because the driver of the m = 4 mode are grid truncation errors which exist over the 
entire domain. 



are dominant and can threaten the integrity of a jet. The higher order "flute" modes (to > 2) corresponding 
to elliptical (to = 2), triangular (to = 3), rectangular (to = 4), etc. modes, do not end up destroying a 
jet, although they can strongly affect the cross-section of the jet and split it into to separate beams [see 
e.g., (Ray 1981)]. Recall as well that the radial structure of modes in jets are of two general types, namely 
surface modes (which are localised towards the surface of the jets) as well as body modes (involving the 
whole body of the jet). Analytical calculations of simple jet models predict that the growth rates of surface 
modes exceeds that of the body modes (see Gill 1965, for discussion). 
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Regardless of what drives the modes, however, once initiated, they will be tracked in 
a physical manner by the numerics, and eventually, the faster-growing m = 1 mode comes 
to dominate the structure. Thus, by t = 200, the jet has responded to the m — 1 mode 
by adopting a helical, or "corkscrew" morphology from z = 30 and beyond. In Fig. 4e, 
the impression is that the collimated jet ends at z — 30, after which the jet has broken up 
into two distinct clumps, one at (y,z) ~ (—7,38), and the other at (y,z) ~ (13,66). In 
reality, the jet remains contiguous throughout the simulation, and the "clumps" are simply 
the intersection of the helical jet with the viewing plane. 

The corkscrew advances in time (i.e., rotates in the same direction as its twist) and by 
t = 240 (Figs. 4f and 5f), the base of the corkscrew has shrunk to z — 15. Three cross sections 
of the corkscrew (and thus more than an entire wavelength) are now contained within the 
primary computational domain [the third "clump" just entering the right-hand boundary 8 
at (y, z) = (—4, 80)], and the centre of the jet has now moved nearly a full jet diameter from 
the disc axis (Fig. 7f). Thus, no part of the jet contains any of the original symmetry axis, 
whence our designation "corkscrew". 

In the cross-sections at z = 30, (Figs. 7-10), one can see the progression from the 
dominance of the quadrantal m = 4 mode (e.g., Figs. 8d, 8e, and 8f) to the helical (m = 1) 
mode (e.g., Figs. 8g and 8h). At this location, much of the effects of the early m = 4 
mode have disappeared by t = 320, while at higher values of z (not shown), the effects of 
the m = 4 mode dissipate significantly earlier (e.g., at z — 50, there is little sign of the 
m = 4 mode by t = 200). In panels f, g, and h of Figs. 7 and 8, steep gradients (indicated by 
"contour bunching") effectively demarcate the cross section of the jet which is predominantly 
elliptical, although evidence of higher order (m > 2) fluting instabilities are apparent in the 
jet profile. Meanwhile, and most significantly, the cross section of the magnetic field (Figs. 
lOg and lOh) shows that a strong axial magnetic field has developed inside the centre (i.e., 
nearer the z-axis) of the displaced jet (c.f., Fig. 7g and 7h), and this acts as a backbone 
providing stability to the jet. 

By t = 320 (Figs. 4g, 5g, and 6g), the corkscrew continues to advance, and has now 
consumed all but the inner several r; of the jet. The wavelength of the corkscrew gradually 
lengthens, but the overall displacement of the jet from the disc axis remains at about one 
jet diameter (e.g., Fig. 7g), with the displacement increasing toward the disc. Between 
t = 240 and t = 320 (e.g., Figs. 7f and 7g), the corkscrew has rotated about 180° in the 



8 Recall that the right-hand boundary is not really a boundary at all, since the "greater computational 
domain" extends to z = 160. Thus, having features and/or material enter from what is ostensibly an outflow 
boundary is quite acceptable, so long as that feature did not originate from the true outflow boundary at 
z = 160. 
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counterclockwise direction. 

By t — 400 (Fig. 7h), the corkscrew has rotated by another 150°, and thus the rotation 
rate is not uniform in time. Further, the rotation rate is more rapid near the base of the 
corkscrew than it is further out. For comparison, at z = 50 (not shown) between t = 240 
and t = 400, the helix advances by only 270°, compared to 330° at z — 30. As a result, 
the pitch of the helix is neither uniform nor constant. Further, severe m > 2 (elliptical and 
higher order) distortions can be seen in the jet profile in panel h of Figs. 7-9, and thus the 
shape of the jet profile is also highly dynamic. Nevertheless, the jet manages to maintain 
collimation throughout the simulation; outflow is never interrupted, nor significantly slowed. 

4.2. How Do Corkscrew Jets Maintain Stability? 

Typical outflow speeds along the jet range from 0.7 to 0.9 in units of the Keplerian 
velocity at r = We argue that if the jet can manage to configure itself in such a way that 
the Alfven speeds are comparable to or higher than unity, then the jet will be sub-Alfvenic, 
and the K-H instabilities will be saturated (thereby satisfying linear stability conditions, 
e.g., Ray 1981). We now show that this Ansatz is precisely satisfied by the behaviour of the 
jet in simulation D. 

Figures 7 and 10 show cross cuts of the density and magnetic field respectively. At all 
later times, t > 240, we see the body of the jet displaced by a jet diameter or so from the 
disc axis. While both the density and normal magnetic field cross sections show well-defined 
peaks, upon careful examination one finds that the peaks are not cospatial (e.g., compare 
Figs. 7g and lOg). This anti-correlation of density and magnetic peaks is found frequently 
in MHD applications, and the present authors have discussed this effect previously (Clarke, 
Norman & Burns 1986; OPII). 

Corresponding to the peak in normal magnetic field are flux loops (shown as vectors in 
Fig. 10), and combined with the normal field, provide the jet with a "backbone" of relatively 
strong helical magnetic field. The density peak is then centrifugally driven toward the outside 
(away from the disc axis) of the backbone as the corkscrew advances in time. 

At every slice along the jet at every epoch after t = 220, the density peaks always cor- 
responds to the location of the highest Alfven Mach number (high density and low magnetic 
field yields a low Alfven speed), and these are typically of order 2 to 2.5. In contrast, the 
centre of the magnetic backbone corresponds to the lowest Alfven Mach number (low density 
and high magnetic field yields a high Alfven speed), and these are typically of order 0.1, or 
less. Averaged over the entire cross section of the jet, the Alfven Mach number is of order, 
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or less than unity, and this is ultimately how the stability of the jet is maintained. 

Thus, our explanation for the maintenance of jet stability goes as follows. At first, when 
the jet is launched, its internal dynamics are dominated by the 2-D symmetry described in 
OPII. The jet has no reason at the very early stages to respond to the m = 1 mode, as it is 
not present until t = 160. Thus, the initial distributions of density, magnetic field strength, 
and Alfven Mach numbers are entirely determined by the nature of the disc and the corona. 

However, the m = 1 mode does grow in time, and the jet must respond. It is not free to 
reduce its velocity, since this is determined by the forces at the base of the jet, and the jet 
begins to buckle under the stresses of the m = 1 mode. But by buckling, the magnetic field 
lines inside the jet are stretched, twisted, and thus strengthened. As this happens, portions 
of the jet become overpressured with magnetic pressure, squeezing out thermal material to 
restore total pressure balance inside the jet. This is ultimately responsible for the separation 
of the density and magnetic peaks, and as a result, portions of the jet attain very low Alfven 
Mach numbers. 

The jet continues to respond to the m = 1 mode by forming a helical structure, and 
as the amplitude of this structure increases, so does the strength of the magnetic backbone 
forming inside the jet, along with the average Alfven speed. At some point, the mean Alfven 
Mach number inside the jet is reduced to or below unity, not because of a reduction in flow 
speed, but because of the increase in Alfven speed. Thus, the jet becomes sub-Alfvenic, and 
the K-H instability is saturated. 

As seen in the simulation, this balance manages to maintain itself with slight variations 
and oscillations for nearly half of the computational time. While we cannot continue the 
present simulation any further than we have because of boundary effects creeping into the 
solution, it seems plausible that this balance should persist indefinitely, giving jets which are 
initially K-H unstable a characteristic helical, or corkscrew morphology. 

Of course, not all jets are observed to be corkscrews, and one immediately asks whether 
there is any way to avoid such a morphology, given the evidence presented in this section 
that corkscrews are a natural configuration for a jet to assume in order to preserve stability. 
Simulation E discussed in §5 offers at least one scenario in which a jet manages to retain 
stability without attaining a large amplitude corkscrew morphology. 
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5. A "WOBBLING" JET 

In this section, we focus on simulation E which we characterise as a "wobbling" jet 
because the centre of the jet never strays more than a jet radius away from the original 
symmetry axis. This simulation uses a smoothed Kepler velocity profile at the inner edge of 
the disk (Fig. 2), in contrast with simulation D which used the cusped profile. Simulation E 
produces a much richer structure in the jet, which prompted us to develop a Fourier trans- 
form analysis package (see http://www.nordita.dk/~ouyed/JETTOOLS) that could reveal 
the quantitative details of the modes that are excited within the jet. 

5.1. Nature of the Outflow 

Figure 11 (and 12) shows isodensity contours of simulation E in the y-z (and x-z) plane 
containing the disc rotation axis at t = 50, 80, 120, 130, 150, 180, 210 and 240. By t = 50 
(Fig. 11a and Fig. 12a), the Alfven wave has moved off the right boundary and the jet 
has propagated nearly half way across the grid. The jet is driven by the field lines that 
are sufficiently displaced radially outwards (i.e., between 1 < r < 5; the jet is hollow), 
which allows the centrifugal acceleration to occur, similar to what was observed in the 2-D 
simulations of OPII. 

The jet at these early times is very stable, and accelerates to maximum velocities of the 
order of 0.7i>K,i- This is somewhat slower than the 2-D counterpart (OPII) where velocities 
of the order of 1.2i>K,i were reached where mass entrainment is not as effective in slowing 
down the jet as it is in 3-D. The flattening (Figs. 11c, lid and lie) and stretching (Figs. 
12c, 12d and 12e) of the jet is evident as well as its response to the kink mode (Figs. 1 If, 
llg, 12f and 12g; see also §5.2.1) before it eventually regains a nearly cylindrical morphology 
centered on the disc axis (Figs, llh and 12h). 

The distortion of the cross-sectional shape of the jet is shown in Fig. 13, which displays 
isodensity contours in the x-y plane located at the vertical line in Fig. 11a. The jet's cross- 
section becomes increasingly elliptical from t — 80 onward, and evidence of higher-order 
fluting modes becomes apparent by t = 120. The highly elliptical cross-section appears 
to break apart into separate streams at t ~ 150. This behaviour is strongly suggestive of 
the non-linear evolution of an m = 2 elliptical mode (see below). We see that this highly 
elliptical, and even bar-like distortion, gradually fades away, so that at t > 200, the jet profile 
appears to be more cylindrically symmetric in the main, with the exception of an obvious, 
one-sided bar-like protrusion that is suggestive of a residual m = 1 helical mode. 

The elliptical cross-section of the jet in Fig. 13 precesses until t = 130, and then appears 
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to remain fixed in position angle until t = 200. This is indicative of equal amplitudes in the 
m = ±2 elliptical modes, discussed further in the next subsection. The precession of the jet 
cross-section resumes after t = 210. 

Figure 14 shows the evolution of 20 magnetic field lines at t — 50, 80, 120, 130, 150, 
180, 210 and 240. The lines were chosen to visualise best the complex dynamics in the jet. 
The two central magnetic field lines (dotted lines) originate from the central compact object 
inside r = r\ and trace the poloidal field lines that ultimately serve as a "backbone" for the 
jet. 

The structure of the jet magnetic field remains well-ordered until t ~ 80. Before this 
time, the inner two field lines remain rather straight, indicating that the axis of the jet is 
quiescent. The inner- most field lines attached to the disk have a clearly helical structure as 
these are associated with the jet itself. The outer-most field lines are beyond the collimated 
outflow and are affected only by the slow rotation of the corona. As such, they are strongly 
poloidal in character. 

The jet's well as its helical structure, become more disorganised at 80 < t < 

210 and execute both a long-wavelength, transverse wandering, as well as shorter-scale, 
disorganised motions. It is clear however, that the jet survives this unstable behaviour and 
appears to resume its initial ordered, regular character at t > 210. Throughout it all, the 
acceleration region of the outflow close to the disk remains largely unaffected. 

Figure 15 shows 20 individual streamlines in the jet at the same epochs as in Fig. 14. 
We note that in general, time-dependent flow, streamlines are not restricted to flow along 
surfaces of constant magnetic flux, and hence do not follow field lines, as is apparent when 
comparing Figs. 14 and 15. Streamlines at t ~ 50, and later at t > 210, show an ordered, 
helical outflow. Between these times however, note the disorder of the flow as various modes 
of instability play themselves out in the evolution of the jet. 

Figure 16 shows the velocity vectors in the same y-z plane as Fig. 11. and further 
illustrates the effect of the kink mode on the body of the jet (Fig. 16e and Fig. 16f) before 
it regains its coherence and cylindrical (one stream) shape in Figs. 16g and 16h. Figure 
17 shows the velocity vectors in the same x-y plane as Fig. 13. One sees that the ordered 
rotation of the jet (a consequence of the fact that it is carrying off the angular momentum 
of the driving disk) is present up to t = 80. This gives way to far more disorganised motion 
between 80 < t < 210. In Figs. 17g and 17h, we see that an ordered, nearly circular rotational 
motion is re-established in this plane. This reflects the behaviour of the streamlines in Fig. 
15. Finally, this evolution of jet rotation is qualitatively similar to what was seen in Fig. 9 
for the corkscrew jet (simulation D, §4). 
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Thus, the jet begins as a stable outflow, destabilises between 80 < t < 210, and then 
resumes some decorum of stability after t = 210. The instability at intermediate times 
appears to be driven by the excitation of a number of discrete K-H modes m > 0. The 
region close to the disk, being characterised by sub-Alfvenic flow, remains reasonably quiet 
throughout the simulation. Thus, as with the corkscrew jet, the wobbling jet manages to 
survive the onset of potentially destructive non-axisymmetric modes, and this is investigated 
in the next sub-section. 



5.2. How do Wobbling Jets Maintain Stability? 

We decompose the radial structure of the jet by performing a Fourier transform (see 
http://www.nordita.dk/~ouyed/JETTOOLS) of the 2-D pressure distributions (not shown, 
but qualitatively similar in appearance to the density distribution shown in Fig. 13). Our 
results are shown in Fig. 18 where we plot the amplitude of a mode with radial wave number 
k r and azimuthal wave mode number, m. (k r is measured in units of (lOr^) -1 ). The grey 
scale ranges from high amplitude (white), to moderate amplitude (grey), and down to low 
amplitude (black). 

5.2.1. Onset of Kelvin- Helmholtz Modes 

As observed for the corkscrew jet, we find that the m = mode is the predominant 
mode at early times, < t < 80, corresponding to when the initial cylindrical symmetry of 
the initial setup survives. Indeed, the high density (and pressure) region of the jet (Fig. 13) 
is nearly circular at these times. The m = mode reappears at later times (t > 210) when 
the jet cross-section is again nearly circular. 

The elliptical, \m\ = 2, modes 9 responsible for the elliptical cross-section of the jet (Fig. 
13), first appear at t = 80, and persist until t ~ 210. . We note that the amplitude of 
both the m = —2 and the m = 2 elliptical modes in Fig. 18, are nearly the same for most 
times, and this freezes the position-angle of the elliptical cross-section in space, as noted 
in §5.1. The \m\ = 2 modes ultimately grow enough to cause the jet to bifurcate between 
120 < t < 180. 



9 There are two senses of rotation for each m-mode of a cylindrical jet. Thus, m > and m < correspond 
to waves that wind around the jet axis in either the same or the opposite sense as the toroidal magnetic field 
respectively. 
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The higher order, \m\ = 4, modes appear slightly after the \m\ = 2 modes starting at 
t = 120, and disappear again at t ~ 180. They manifest themselves by giving the cross- 
section a marked rectangular appearance, such as in Figs. 13c and 13d. As before, these 
are excited, in part, by the initially cylindrical symmetric atmosphere rotating within a 
quadrantally symmetric grid. 

The \m\ — 1 modes makes their first strong debut (relative to the amplitude of the flute 
and pinch modes) rather late in the evolution of the jet, at t ~ 150. The helical mode is of 
comparable amplitude to the elliptical modes, so the jet's cross-section still remains rather 
elliptical in shape. 

It is clear from Fig. 18 that for 150 < t < 210 the jet's evolution is dominated by the 
m = l and m — — 1 modes, with the higher order flute modes diminishing in importance. We 
also note that the amplitude of the m = 1 mode is always greater than that of its m — —1 
counterpart. This is why the late- appearing, bar-like protrusion noted earlier undergoes its 
slow rotation. 

A global view of Fig. 18 shows that the various modes that are strongly in play at 
t = 120, gradually damp out in their amplitude. By the end of our simulation, at t — 240, 
there is a bit of activity in the m — 1 mode. This "spectrographic" image of the jet's 
evolution shows that the jet has survived the instability, and that the potentially unstable 
modes all damp down with time. The jet has all but evaded the most threatening instability. 
A second general trend apparent in Fig. 18 is that the radial wavenumber of the unstable 
modes becomes smaller, that is k T — > (e.g., compare Fig. 18b with 18h). The increasing 
radial wavelength of the modes with time shows that the jet is regaining its coherence from 
a beam that was broken into two separate streams to a single coherent stream again at the 
end of the simulation. 



5.2.2. Stabilising the Wobbling Jet: Transition to Sub-Alfvenic Flow 

Figure 19 shows the time evolution of the Alfven Mach number (Ma) along the innermost 
magnetic field line (r = n) at 50 < t < 240. In these panels, the value of Ma is plotted as a 
function of the position s along the field line. The vertical dashed line in any frame indicates 
the point sa along the field line at which the flow reaches the Alfven point (where Ma = 1). 

Up to t — 80, the maximum value of Ma continues to increase beyond the Alfven point, 
to a maximum of Ma = 4. The position of the Alfven point is at s a ~ 12-15 at these earlier 
times. At t > 80, we see that the maximum value of Ma decreases systematically with time. 
At 150 < t < 210, the velocity of the flow on this entire field line is sub-Alfvenic. We note 
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that t = 150 in Fig. 19 is approximately the time where the \m\ = 2 modes begin to fade 
out. 

This decrease in M A along the field line begins at almost the same time as the appearance 
of the non-axisymmetric modes in the jet. We note that the reappearance of super-Alfvenic 
flow at t > 210 is the time in which only low level \m\ — 1 modes are left, as discussed 
above. Our results show therefore that the appearance of unstable modes in the jet is 
directly correlated with the Alfvenic Mach number of the jet. 

The instability breaks out when the jet becomes super-Alfvenic. However, the jet Alfven 
Mach number is systematically reduced as a consequence of the unstable modes. The modes 
begin to lessen in strength as soon as the Mach number decreases below unity. This is 
excellent evidence that flow stabilisation occurs as a consequence of the mechanism that 
de-stabilises the jet in the first place. Once again, we see that the jet is self-regulating in 
the sense of being able to adjust its structure to preserve conditions of stability. 

Figure 20a shows the time evolution of the mean average toroidal magnetic field in the 
jet. The toroidal field builds in amplitude while the jet is stable, reaching a peak at t — 50. It 
decreases with time beyond this until t ~ 130 and remains at this lowest value until t ~ 180 
whereupon it starts a steady and monotonic rise to a maximum value achieved at the end of 
the simulation. 

The overall energetics of the jet are shown in Fig. 20b. Despite the fact that decreases 
between 50 < t < 130, the total magnetic energy systematically increases because of an 
increase of the strength of the poloidal field component caused by the stretching of field lines 
by the unstable \m\ — 1 and higher order fluting modes. 

The fact that the magnetic energy taps into the bulk kinetic energy of the flow (which 
includes energy stored in the K-H modes) is apparent in Fig. 20b, where between t = 130 
and t = 170, E m rises as E k falls. Just as in the case of the corkscrew jet, the sheared 
velocity field in the jet is the ultimate source of the energy tapped by the unstable modes 
and transferred into poloidal field, stabilising the jet against the \m\ = 1 modes. 

After t = 180, the kinetic energy grows at the expense of the magnetic tension in 
the poloidal field, and the jet is efficiently accelerated to super-Alfvenic flow once again. 
This quasi-periodic transfer of energy between kinetic and magnetic fields suggests that 
the jet may oscillate between quasi-stable and quasi-unstable epochs as it propagates into 
its surrounding environment. Unfortunately, our simulation could not be carried out long 
enough to confirm this hypothesis, because effects of the outflow boundary conditions were 
beginning to creep into the primary computational domain. 
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6. DISCUSSION AND CONCLUSIONS 

We have presented results from a numerical study of a variety of different simulations 
of 3-D winds from accretion disks threaded by vertical field lines. Although, our present 
simulation is far from the high resolution used in previous 2-D simulations (OPI and OPII), 
we find that the acceleration phase of the jet is very similar. 

Our set of simulations suggests a general stabilisation mechanism independent of how 
the jet is perturbed, and on the nature of the boundary conditions employed. The most 
fascinating result of our simulations is that jets in 3-D, while manifesting the expected 
unstable modes long discussed in the literature, remain stable. This remarkable behaviour 
can be traced to the self-limiting action of the instability itself. The jet becomes unstable 
in regions of high Alfven Mach number, as both analytical and numerical studies of much 
simpler jet systems have suggested (e.g., Ray 1981; Hardee & Rosen 1999). The appearance 
of the \m\ = 1 modes pump energy into the poloidal magnetic field, causing the jet Alfven 
Mach number to fall below unity and stabilise the jet. The unstable modes die away, and 
the flow can once again start to accelerate to Alfven Mach numbers greater than unity. 

The differences between simulations D and E can be traced to the different profiles 
imposed at the accretion disc (Fig. 2). In short, the cusped profile used in simulation D moves 
more specific angular momentum per unit time onto the grid than the smoothed profile of 
simulation E. Thus, the m = 1 mode is highly dominant in simulation D, driven by the 
high angular momentum of the inner parts of the corona. In contrast, this driver is subdued 
in simulation E, and other modes of the K-H instability are manifest. The simple m — 1 
dominance in simulation D leads to the corkscrew morphology in which the jet achieves a 
balance between the centrifugal and magnetic stresses (very similar to OPII explanation for 
the m=0 mode in the 2-D jet, and the knot generator). With lots of modes playing a role, 
the dynamics are much more complex in simulation E. 

Above all, the critical point is that despite the differences in boundary conditions, both 
jets manage to strike a balance and avoid disruption by the many modes dumping energy 
into the magnetic field. 

The development of a corkscrew morphology is interesting in its own right. While we 
are quick to emphasise that the scale of our simulation corresponds only to the very closest 
regions to the compact object, it is nevertheless noteworthy that the two closest examples of 
extragalactic jets, Centaurus A (Clarke, Burns, & Feigelson 1986) and M87 (Biretta, Zhou 
& Owen 1995) show definite side-to-side oscillations and/or helical morphologies in their 
structures. This behaviour is even clearer in the observations of wiggling optical jets from 
YSOs such as HH47 (Heathcote et al. 1996) and GGD 34 (Gomez de Castro & Robles 1999). 
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These jets are observed in optical forbidden lines, typically SII. As noted by Heathcote et al. 
(1996), supersonic jets tend to move ballistically making bending difficult to achieve. We note 
however that the relevant quantity that constrains the behaviour of magnetised flows is their 
Alfven Mach number, which our simulations show adjust to near unity. Our simulations occur 
on scales that are about a decade smaller than the HST can resolve for nearby protostars. 
Nonetheless, the fact that some of our jet simulations settle into a wobbling or corkscrew 
configuration probably says something interesting about the larger scale behaviour as well. 

Finally, we have seen that our simulated jets regain a reasonably stable configuration at 
late times. Indeed, various physical quantities in the jet, such as the strength of the toroidal 
field, the jet Mach number, etc., are rather close to their values at earlier times, before the 
unstable flute modes make their appearance (e.g., at t < 50). This strongly suggests that 
the jet would ultimately undergo another burst of unstable behaviour after it accelerates 
material to sufficiently high Alfvenic Mach number. Testing this conjecture requires a much 
larger simulation wherein one doubles or triples the size of the box to ensure that most of 
the body of the jet stays in the computational domain. It seems clear however, that the jet 
will undergo this kind of episodic behaviour, which in turns has interesting consequences for 
the appearance of multiple bow-shocks and other episodic phenomena that are seen in real 
jets. 
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A. 



BOUNDARY CONDITIONS 



A.l. Truncating the Keplerian profile 



A Keplerian profile, 



v K (r) 



1 



(A.l) 



y/r 



is impractical to implement over all r on a numerical, Cartesian grid. In particular, the 
singularity at r = must be avoided, and at the outer limits of the grid, due attention must 
be paid to the rotation in the grid "corners" . 

The singularity can be handled most conveniently by setting ity(r) = for r < r^. This 
solution, depicted in Fig. 2, results in a cusp at r = r^, which may introduce numerical 
concerns of its own. Alternatively, one may introduce a transition region over which the 
Keplerian profile is continuously and smoothly connected to the v^{r) = solution inside r^, 
also shown in Fig. 2. Both options were used in the simulations discussed in this paper, and 
the mathematical form of the smoothing function used is given below [equation (A. 2)]. 

As for the outer regions of the grid, consider the corner of the computational domain in 
the (+x, +y) quadrant, as depicted in Fig. 21. Here, two outflow boundaries meet, namely 
the +x-boundary and the +y-boundary. If rotation about the z-axis is positive (in the 
sense of the right-hand rule), then near the (+x,+y) corner, the sense of rotation causes 
material to flow out across the +y-boundary, which as an outflow boundary it can handle. 
On the same token, material should flow in through the +x-boundary, which as an outflow 
boundary, cannot happen. Thus, a vacuum is created in the corners, adversely affecting the 
rest of the computational domain within a sound-crossing time. 

For many reasons, it is impractical to suggest that the +x boundary be redesignated 
as an inflow boundary. This would break the symmetry that must exist between the +x- 
and +?/-boundaries, and would require inflow values to be specified a priori all along the 
+a;-boundary, which is not possible without the full time- dependent solution to the problem! 

This problem doesn't exist in cylindrical coordinates, of course, since there are no cor- 
ners. We therefore try to mimic this desirable property of cylindrical coordinates by trun- 
cating the ity(r) profile at an "outer radius" r , where r is the radius of a cylinder which 
completely contains the primary computational domain, yet lies well inside the greater com- 
putational domain (Table 1, and Fig. 21). Thus, if there is no rotation in the corners, the 
problem of evacuating the corner regions is avoided. 

If the truncation at r D is sudden, then the outer boundary of the rotating cylinder will 
be ragged (being resolved on a Cartesian grid). This will create a large and undesirable 
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numerical viscosity along the outer surface of the cylinder, which can be largely avoided if 
we use a smooth truncation instead. Thus, we impose a continuous and smooth profile for 
v^{r) between r Q and r max , where r max > r Q but still less than the maximum value along the 
x- and y-axes (for simulations A, B, C, and E, we used r Q = 22 and r max = 28, while for 
simulation D, we used r Q = 30 and r max = 38). 

Therefore, the form of the azimuthal velocity profile, v^r), used in these simulations 
which truncates the Keplerian profile [equation (A.l)] at both the inner and outer regions of 
the grid is given by: 

r < n, 
r 1 <r < r 2 , 

r 2 <r < r , (A. 2) 

To — T *^ ''maxi 
T — '"max ? 

where all variables are scaled as described in §2.1. For simulations A-D, r\ = r 2 = r i; and 
there is no region in which the sin 2 smoothing function is applied. For simulation E, r\ = 
and r 2 = 2r\ and the sin 2 solution raises the profile smoothly from to the Keplerian value. 
For all simulations, the cos 2 function allows the profile to drop off smoothly between r Q and 
''max, arriving at the minimum value (0) with zero slope. 

Of course, once flow has advanced onto the grid by a few zones, the dynamics of the flow 
will have altered the velocity field significantly from whatever was specified at the boundary. 
Thus, which functions are used to smooth the profiles should not be critical. 
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A. 2. Breaking the quadrantal symmetry 

When attempting to break a geometric symmetry, experience has shown that it is prefer- 
able to impose perturbations to the velocity fields, rather than directly to the pressure or 
density distributions. Pressure perturbations send sound waves in all directions, which in a 
highly dynamic situation, can steepen into shocks and affect outlying regions more than the 
original perturbations may have intended. 

Thus, we break the imposed quadrantal symmetry by introducing a non-azimuthal com- 
ponent to the velocity profile of the disc. We do this in two ways. The first, and perhaps 
most obvious way, is to offset the point about which the velocity profile in (A. 2) is evaluated 
from the centre of the gravitational potential well. Then at some prescribed time, tj, later, 
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the origin for the velocity profile is "jerked" to another location. This is implemented by 
replacing r = \J x 2 + y 2 in equation (A. 2) with £, where £ is given by: 

e = (a.3) 

[ - (y - t > tj, 

where 5r is a few or a few tens of percent of (it does not seem to matter), and where 
we used tj = 10 (in scaled time units described in §2.1), though again, the precise time the 
origin is "jerked" does not matter. Simulations B and E were "jerked" in this manner to 
break the quadrantal symmetry. 

For simulations C and D, the disc was "wobbled" by introducing a time-dependent radial 
component to the disc velocity profile. Thus, we use a disc velocity of the form: 

v = v r (r, <j),t)r + v < p(r)4> (A.4) 

where v<f,(r) is given by equation (A. 2), and where the radial component, v r (r,<f>,t), is the 
perturbation and is given by 

v r (r, <f>,t) — - cos(0 + out). (A. 5) 

As usual, the azimuthal coordinate, 0, is taken counter-clockwise from the +x-axis, where 
the z-axis is the disc axis. We used e = 0.02 (a 2% perturbation at r = Ti = 1), and 
uj — 0.1 (one tenth the Keplerian angular velocity at r = Tj = 1). The introduction of 
the radial velocity component distorts the otherwise circular velocity profile imposed by 
Vj,{r) alone into a precessing elliptical profile, with the origin of the gravitational potential 
located at the centre of the ellipse. The precession is introduced by virtue of the explicit 
time-dependence in (A. 5), without which, the velocity profile would still possess quadrantal 
(though not azimuthal) symmetry. 

We close this section by noting that "jerking" and "wobbling" the disc seems to give 
qualitatively identical results. Thus, how one chooses to break the quadrantal symmetry is 
a matter of taste, so long as the symmetry is truly broken. 



A.3. Enforcing V ■ v = at the disc 

Another consequence of using Cartesian coordinates is the numerical introduction of a 
divergence in what analytically is a solenoidal vector field. Physically, the velocity profile 
given by equation (A.4) is solenoidal at any given instant of time. This may be verified 
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by noting that the velocity streamlines form closed loops, or by evaluating the divergence 
directly. 

However, as the disc is resolved on a Cartesian grid, the x- and ^/-components of v must 
be evaluated, and if this is done by direct coordinate transformation (e.g., v x = v r cos(p — 
v<psm(p, etc.), numerical truncation errors can generate significant non-zero divergences in a 
time short compared to the duration of the simulation. 

First, we establish the need to preserve the solenoidal nature of the velocity profile to 
within machine roundoff errors, and then we show how this can be done. Consider the task of 
preserving B z to a constant value on the surface of the disc, as required by these simulations. 
From the induction equation, we have 

dB z = de ]L _de x _ 
dt dx dy ' 

with e x and e y being the x- and y-components of the so-called e.m.f. (Evans & Hawley 1988), 

— * 

defined as e = v x B. Thus, for v z = 0, 

e x = v y B z - v z B y = v y B z , (A. 7) 

e y = v z B x - v x B z = -v x B z (A. 8) 

For B z to stay constant on the disc surface, equations (A. 6), (A. 7), and (A. 8) require that 



dB z 



dt 



z=0 



'dv x dv y 

—£>z I "5 r tt- 

ox ay 



= -B Z W ■ v 

2=0 



2=0 



0, (A.9) 



Thus, if V • v 7^ to machine round-off errors (as would be the case if one naively evaluates 
v x and v y from the components of v), B z will evolve in time with profound and disastrous 
consequences. In our case, we noted that the growth of B z was most severe near the inner 
boundary of the disc, causing the inner region of the atmosphere to be evacuated via nu- 
merically driven "jets". With the Alfven speed (B/y/p) unbounded, the Alfven time step 
vanished and the simulation ground to a halt. 

The fix is simple, and recognisable by anyone who has faced the task of setting up a 
numerically solenoidal magnetic field. Consider the vector quantity q given by v = V x q. 
Note that q is to the velocity field what the vector potential is to the magnetic field. In 
cylindrical coordinates, if v = v r r + v <p 4> [e.g., equation (A. 4)], then q = q z z, and; 

1% dq z 

v * = rW v * = -~dV (A - 10) 
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and 



For the purposes of finding q, let us assume that q z is separable. Thus, let q z {r,<f>) = 
R(r) + $(0). Then, inverting equations (A. 10), we get 

nr 

R(r) = - v^{r')dr' (A.ll) 

$(0) = [* v r (<f>') d(j>' (A. 12) 

Jo 

where is given by equation (A. 2) and v r is given by equation (A. 5). Some of the cases 
in equation (A. 2) are integrable analytically, and some require numerical integration. Re- 
gardless, for any given value of r, one can find R(r) from equation (A.ll). Meanwhile, v r 
integrates trivially, giving 

q z = R(r) + e [sin(0 + ut) - sin(cjt)] (A. 13) 

Now, the z-components of vectors are invariant under the transformation between cylin- 
drical and Cartesian coordinates. Thus, we may consider equation (A. 13) to give the z- 
component of q in Cartesian coordinates, and evaluate the x- and y-components of the 
truncated, perturbed Keplerian velocity profile from q z . Thus; 

v *=v Vy = ~^ (A - 14) 

In particular, for a staggered mesh such as that employed by the ZEUS family of codes 
(Clarke, 1996), q z as a quantity would be centred on the zone-edges parallel to the z-axis 
(Fig. 22). Thus, a y-difference of q z centres v x on the x-face (as required on the staggered 
mesh) while an x-difference of q z centres v y on the y-face (again, as required on the staggered 
mesh). Thus, the difference form of equations (A. 14) is 

Vx(ij) = riU+iHM , (A.15) 

vy{ij) = _ ,,(.• + (A 16) 

where we now write the variables as functions of their coordinate indices, rather than 

of the coordinates, (x,y), themselves. 

With the numerical divergence given by 
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it is left as a straight-forward exercise to show that the numerical divergence of v will be zero 
to within machine round-off errors. Further, when averaged to the zone centres and combined 
quadratically, the velocity component in equations (A. 15) and (A. 16) reproduce the original 
velocity profile (equation A. 4) to within a few percent, depending on the coarseness of the 
grid chosen. 

A. 4. Magnetic Outflow Conditions 

To be frank, it is still unclear what the "optimal" magnetic outflow boundary conditions 
should be. In the ZEUS-2D simulations described in OPI and OPII, for example, both 
(the only component of the e.m.f. tracked in ZEUS-2D, and used to update the poloidal 
components of the magnetic field), and are kept constant across outflow boundaries. 
Note that setting e<j, constant across the boundary is not at all the same as setting the 
poloidal field components themselves constant across the same boundary. For example, if 
B r is initialised everywhere to zero and is set constant across outflow boundaries (as in 
OPII), B r will forever remain zero outside a ^-outflow boundary (because de^/dz remains 
zero across the boundary) even if B r just inside the boundary were to become non-zero. On 
the other hand, if B r were set constant across the outflow boundary, the value of B r outside 
the z-boundary would change with the values just inside the grid. 

Thus, the toroidal and poloidal magnetic field components are treated differently across 
an outflow boundary in ZEUS-2D, and yet there were no physical principles used to justify 
this (Jim Stone, private communication). Still, when one tries the obvious alternative of 
setting the poloidal components of the magnetic field constant across an outflow boundary, 
severe numerical errors occur at the boundary yielding an unphysical build-up of magnetic 
stresses which evacuate the boundary zones. The combination of large magnetic field and 
low density results in very high Alfven speeds and thus vanishingly small Alfven time steps, 
grinding the simulation to a halt. 

Unfortunately, there is no obvious way to extend the ad hoc prescription for magnetic 
outflow boundary conditions used in 2-D cylindrical coordinates (namely, maintaining con- 
tinuous B<j, and across boundaries) to 3-D Cartesian coordinates. For the present, we 
set all magnetic field components constant across outflow boundaries, and then impose floor 
values on the density near the outflow boundary to prevent the Alfven time step from van- 
ishing. To avoid the bad densities and ensuing dynamics from corrupting the solution as 
soon as the outflow reaches the boundary, we have pushed the actual outflow boundary well 
away from the primary region of interest, and filled the intervening regions with increasingly 
coarse zones (Fig. 1, §3.2). Eventually, the effects of the magnetic stresses building up on the 
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distant outflow boundaries make their presence felt on the primary computational domain, 
and it is at this point that the simulation is stopped. While not an ideal solution, it is a 
practical one, without which, the simulations presented herein could not have been evolved 
as far as they were. 
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Fig. 1. — A schematic of the y-z slice of the 3-D grid for simulations A, B, C, 
and E that contains the disc (z) axis. The x-axis, which is zoned identically to 
the y-axis, points into the page at z — y — 0. Inflow conditions are imposed 
at the left boundary (y-axis) which includes the accretion disc between r» and 
r max (designated by the heavy black line). All other boundaries (dashed) are 
designated as "outflow" boundaries. In this schematic, only every other zone is 
represented. 
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Fig. 2. — Profiles for v^r) in the disc [given by equation (A. 2) in Appendix A] for 
simulations A-C (cusp at r = r^) and E (smoothed in the region < r < 2rj). 
The f</,(r) profile used in simulation D is the same as the cusped profile, except 
it extends to r = 40rj, with the cutoff occurring between 30rj < r < 38?v The 
dramatic differences between simulations D and E can be largely attributed to 
the different profiles for v$(r). 
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a) 



b) 




Fig. 3. — False colour representation of a) J pdl and b) J V • vdl for simulation 
D. The disc is on the left side (not visible), and flow is generally from left to 
right. Colours are arranged spectrally from blue to red to represent low and high 
values of the variable. 
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Fig. 4. — 2-D contour slices of the density on the y-z plane [where the z-axis 
(horizontal) is the disc axis] for simulation D shown at t — a) 20, b) 60, c) 100, 
d) 160, e) 200, f) 240, g) 320, and h) 400. The +x-axis, located at y = z = 0, 
points into the page. H and L indicate local maxima and minima respectively. 
The vertical line in panel a) (at z = 30) indicates the location of the cross- 
sectional slices in Figs. 7-10. The y-axis extends to ±20. 
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Fig. 5. — 2-D contour slices of the (toroidal) magnetic field component going into 
(dashed contours) or coming out of (solid contours) the page on the y-z plane 
[where the z-axis (horizontal) is the disc axis] for simulation D shown at the same 
times as Fig. 4. H and L indicate local maxima and minima, respectively. The 
y-axis extends to ±20. 
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Fig. 6. — 2-D vector slices of the poloidal velocity components [within the y-z 
plane where the z-axis (horizontal) is the disc axis] for simulation D shown at 
the same times as Fig. 4. The y-axis extends to ±20. 



-42- 



a) t = 20 b) t = 60 c) t= 100 d) t = 160 









/ i i i | i i i \ 








^/ 1 1 | 1 1 


?» 0.0 


^^^^ 




\ i i i~i i i i / 




^\^i~^~T~ — 1 — \ /\ 




^^^n. i i i i i 



0.0 0.0 0.0 0.0 



x x x x 



e) t = 200 f)t = 240 g) t = 320 h) t = 400 




0.0 0.0 0.0 0.0 



x x x x 



Fig. 7. — 2-D contour slices of density on the x-y plane at z = 30 (corresponding 
to the vertical line in Fig. 4a) for simulation D shown at the same times as Fig. 4. 
H and L indicate local maxima and minima respectively. The x and y axes extend 
to ±20 (i.e., the entire primary computational domain). From this vantage, the 
z-axis points out of the page from the centre of each plot, and the sense of disc 
rotation is counter-clockwise. At early times, the circular contours reflect the 
initial spherically-symmetric atmosphere, but by t = 200 (panel e), the m = 4 
mode arising from the quadrantal symmetry of the grid is apparent. At higher 
times, the m = 4 mode gives way to the m = 1 mode, forcing the jet axis well 
off the z-axis of the grid. 
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Fig. 8. — 2-D contour slices of the Alfven Mach number (Ma) on the x-y plane 
at z = 30 for simulation D shown at the same times as Fig. 4. Particularly from 
t = 200 and on, the, maxima of Ma near the core of the jet correspond to the 
density maxima seen in Fig. 7. 
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Fig. 9. — 2-D contour slices of the normal velocity (v z ) with poloidal velocity 
vectors {v x x + v y y) superimposed on the x-y plane at z = 30 for simulation D 
shown at the same times as Fig. 4. Dashed contour lines indicate flow into the 
page. Maximum normal velocities correspond to maxima in both density (Fig. 
7) and normal magnetic field (Fig. 10). 
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Fig. 10 . — 2-D contour slices of the normal magnetic field (B z ) with poloidal 
magnetic field vectors (B x x + B y y) superimposed on the x-y plane at z = 30 for 
simulation D shown at the same times as Fig. 4. At t = 20, the toroidal magnetic 
field wraps clockwise about the z-axis, as expected from the torsional Alfven wave 
launched by the counter-clockwise rotation of the disc. At later times, a magnetic 
"spine" (compact circular feature in panels f, g, and h) develops. The density 
maxima (Figs. 7g and 7h) are offset from the spine away from the z-axis, as 
though the centre of mass of the jet is being driven centrifugally outwards from 
the magnetic spine by the rotation of the jet. 
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Fig. 11. — 2-D contour slices of density on the y-z plane [where the z-axis (hori- 
zontal) is the disc axis] for simulation E shown at t — a) 50, b) 80, c) 120, d) 130, 
e) 150, f) 180, g) 210, and h) 240. The +x-axis (located at y = z = 0) points 
into the page. The vertical line in panel a) (at z = 25.0) indicates the location 
of the cross-sectional slices in Figs. 13 and 17, and used to create Fig. 18. 
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Fig. 12. — 2-D contour slices of density on the x-z plane [where the z- axis (hori- 
zontal) is the disc axis] for simulation E shown at t — a) 50, b) 80, c) 120, d) 130, 
e) 150, f) 180, g) 210, and h) 240. The +y-axis (located at x = z = 0) points 
into the page. 
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Fig. 13. — 2-D contour slices of density on the x-y plane at z — 25.0 (correspond- 
ing to the vertical line in Fig. 11a) for simulation E shown at the same times as 
Fig. 11. 
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Fig. 14. — Snapshots of 20 magnetic field lines for simulation E shown at the same 
times as Fig. 11. The two central magnetic field lines (dotted lines) originate on 
the central compact object (illustrated by the semi-sphere to the left). They are 
not attached to the disk surface (r < r^) and do not rotate. The disc axis is along 
the diagonal of the frame (on a 45° angle). 
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Fig. 15. — Snapshots of 20 streamlines for simulation E shown at the same times 
as Fig. 11 and with the same orientation as Fig. 14. The density isosurface is 
shown to illustrate the regions of high pressure. An appropriate density isosurface 
was chosen which best shows the section of the disk which participates in the 
outflow (r < 5rj). 
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Fig. 16. — 2-D vector plots of poloidal velocity in the y-z plane [where the z-axis 
(horizontal) is the disc axis] for simulation E, shown at the same times as Fig. 
11. Only the inner 2/3 of the plane [(y,z) = (—10 : 10,0 : 60)] is shown. The 
maximum vector length is 0.5Vk,%- 
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Fig. 17. — 2-D vector plots of poloidal velocity in the x-y plane at z — 25.0 for 
simulation E, shown at the same times as Fig. 11. Only the inner half of the 
plane [(x,y) = (—10 : 10,-10 : 10)] is shown. The maximum vector length is 
0.5V^ 




Fig. 18. — Amplitudes of modes with radial wave number (k T ) and azimuthal 
wave number (m) for the x-y pressure slice at z — 25.0 for simulation E shown at 
the same times as Fig. 11. The grey scale ranges from white (high amplitudes) 
to black (low amplitudes). 
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Fig. 19. — Alfven Mach number (Ma) along the innermost magnetic field line of 
simulation E shown at the same times as Fig. 11. The s-axis is the distance along 
the field line (arc length). The vertical dashed lines mark the location, sa, of the 
Alfven point. 
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Fig. 20. — a) The time evolution of the mean toroidal magnetic field integrated 
over the primary computational domain of simulation E. b) The time evolution 
of the bulk magnetic (E m ), kinetic (E k ), and thermal (P) energies, integrated 
over the primary computational domain of simulation E. 
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Fig. 21. — An x-y slice of the grid (for simulations A, B, C, and E) as viewed from 
above the +z-axis. Rotation of the disc in the counter-clockwise direction would 
require material to flow in across the +x-boundary, and then out again across the 
+!/-boundary. Since both boundaries are outflow, the first requirement is impos- 
sible. Thus, the Keplerian profile of the disc is reduced to zero smoothly between 
r Q and r max [equation (A. 2)], preventing the corners from being evacuated. As in 
Fig. 1, only every second zone is indicated. 
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Fig. 22. — In the staggered mesh, variables are not all cospatial. In particular, 
scalars such as the density are located at the zone-centres, while primary vec- 
tor components such as the velocity are face-centred as shown. Derived vector 
components, such as the "velocity vector potential" (q z ) are edge-centred. 



